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Abstract 


In  this  paper  we  employ  the  periodic  unfolding  method  for  simulating  the  electromagnetic  held  in  a 
composite  material  exhibiting  heterogeneous  microstructures  which  are  described  by  spatially  periodic  pa¬ 
rameters.  We  consider  cell  problems  to  calculate  the  effective  parameters  for  a  Debye  dielectric  medium  in 
the  cases  of  circular  and  square  microstructures  in  two  dimensions.  We  assume  that  the  composite  materials 
are  quasi-static  in  nature,  i.e.,  the  wavelength  of  the  electromagnetic  held  is  much  larger  than  the  relevant 
dimensions  of  the  microstructure. 
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1  Introduction 


In  this  article  we  study  the  behaviour  of  the  electromagnetic  field  in  a  material  presenting  hetero¬ 
geneous  microstructures  (composite  materials)  which  are  described  by  spatially  periodic  param¬ 
eters.  We  will  subject  such  composite  materials  to  electromagnetic  fields  generated  by  currents 
of  varying  frequencies.  When  the  period  of  the  structure  is  small  compared  to  the  wavelength, 
the  coefficients  in  Maxwell’s  equations  oscillate  rapidly.  These  oscillating  coefficients  are  difficult 
to  treat  numerically  in  simulations.  Homogenization  is  a  process  in  which  the  composite  mate¬ 
rial  having  a  microscopic  structure  is  replaced  with  an  equivalent  material  having  macroscopic, 
homogeneous  properties.  In  this  process  of  homogenization  the  rapidly  oscillating  coefficients  are 
replaced  with  new  effective  constant  coefficients.  The  primary  objective  of  homogenization,  i.e. ,  of 
the  micro-macro  approach,  is  to  replace  a  system  with  periodically  varying  coefficients  by  a  limiting 
homogeneous  system  that  facilitates  computation.  The  approach  that  we  take  here  is  based  on  the 
periodic  unfolding  method  presented  in  [8,  9,  7].  We  first  mention  other  efforts  on  homogenization 
of  Maxwell’s  equations. 

In  [25],  a  method  based  on  spectral  expansions  for  Maxwell’s  equations  is  presented,  which 
utilizes  eigenvectors  of  the  curl  operators  combined  with  the  microscopic  description  of  the  material. 
The  homogenized  material  is  represented  using  mean  values  of  only  a  few  eigenvectors.  This  method 
relies  on  the  material  being  lossless,  in  which  case  Maxwell’s  equations  can  be  associated  with  a 
self-adjoint  partial  differential  operator.  However,  most  materials  usually  have  losses  due  to  a 
small  conductivity  or  dispersive  effects,  which  renders  the  corresponding  operator  in  Maxwell’s 
equations  non-selfadjoint.  In  [24]  the  authors  use  a  singular  value  decomposition  for  analyzing 
non-selfadjoint  operators  that  arise  in  Maxwell’s  equations.  They  expand  the  electromagnetic  field 
in  the  modes  corresponding  to  the  singular  values,  and  show  that  only  the  smallest  singular  values 
make  a  significant  contribution  to  the  total  field  when  the  scale  is  small.  Using  this  approach  they 
find  effective,  or  homogenized,  material  parameters  for  Maxwell’s  equations  when  the  microscopic 
scale  becomes  small  compared  to  the  scale  induced  by  the  frequencies  of  the  imposed  currents. 
In  [12],  the  authors  compare  two  different  homogenization  methods  for  Maxwell’s  equations  in 
two  and  three  dimensions.  The  first  method  is  the  classical  way  of  determining  the  homogenized 
coefficients  [9],  which  consists  of  solving  an  elliptic  problem  in  a  unit  cell.  The  second  method 
based  on  spectral  expansions  is  described  in  [25].  In  [16],  the  author  presents  an  overview  of 
the  homogenization  of  anisotropic  materials  at  fixed  frequency  using  the  concept  of  two-scaled 
convergence  [20,  1].  The  homogenized  electric  and  magnetic  parameters,  the  relative  permittivity 
and  the  relative  permeability,  respectively,  are  found  by  suitable  averages  of  the  solution  of  a  local 
problem  in  the  unit  cell.  In  [13]  a  homogenization  technique  for  harmonic  Maxwell  equations  in 
a  composite  periodic  medium  is  presented.  See  also  [26,  17,  21]  for  some  other  homogenization 
methods  for  Maxwell’s  equations. 

In  this  paper  we  use  the  periodic  unfolding  method  introduced  in  [8]  in  the  abstract  framework 
of  stationary  elliptic  equations.  The  study  in  this  paper  considers  constitutive  laws  that  take 
into  account  bianisotropy,  chiral  symmetry,  thermal  and  memory  effects.  The  homogenization 
procedure  yields  a  limit  constitutive  law  different  from  the  original  one  wherein  the  convolution 
operator  that  accounts  for  memory  effects  is  replaced  by  a  more  complex  Hilbert-Schmidt  operator. 
We  refer  the  reader  to  [6,  5]  for  the  relevant  theory.  In  the  following  sections,  we  present  the 
electromagnetic  problem  that  is  of  interest  to  us  and  then  set  up  the  corresponding  homogenized 
problem  to  be  solved.  A  comparison  is  made  between  the  effective  parameters  obtained  by  the 
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exact  homogenization  method  presented  here  and  those  computed  by  traditional  mixture  formulae, 
such  as  the  Maxwell  Garnett  formula,  the  Bottcher  mixture  rule  or  Bruggeman  formula,  some  of 
which  are  based  on  physical  arguments  [22]. 

Our  efforts  here  are  motivated  by  use  of  electromagnetic  interrogating  signals  (possibly  in  the 
Terahertz  range)  for  detection  of  defects  [4,  14]  in  the  insulating  foam  on  the  fuel  tanks  of  the 
NASA  space  shuttles.  Defects  in  the  foam  are  believed  to  contribute  to  the  problem  of  separation 
of  the  foam  during  liftoff,  resulting  in  significant  damage  to  and  possibly  subsequent  destruction  of 
the  space  vehicle  itself.  The  low  density,  closed  cell  foam  is  a  very  complex  heterogeneous  material 
[19].  It  is  a  polyurethane-type  foam  composed  of  five  primary  substances:  polymeric  isocyanate,  a 
flame  retardant,  a  surfactant,  a  blowing  agent,  and  a  catalyst.  The  surfactant  controls  the  surface 
tension  of  a  liquid  and  thus  cell  formation.  The  blowing  agent  creates  the  foam’s  cellular  structure 
by  creating  millions  of  tiny  bubbles  or  foam  cells.  As  a  first  approximation,  we  consider  here 
materials  with  periodic  gas  filled  cells  surrounded  by  a  matrix  of  polyurethane  type  nonmagnetic 
material.  The  dielectric  properties  (permittivity,  conductivity,  etc.)  vary  substantially  between 
the  cellular  and  matrix  materials,  leading  to  highly  oscillating  coefficients  in  the  Maxwell  system 
describing  propagation,  reflection  and  dispersion  of  the  electromagnetic  fields  resulting  from  the 
interrogating  probes. 


2  Maxwell’s  Equations  in  a  Continuous  Medium 


We  employ  Maxwell’s  equations  for  a  linear  and  isotropic  medium  in  a  form  that  includes  terms 
for  the  electric  polarization  given  by 

3D 

(i)  —  +Jc-VxH  =  J,  in  (0,T)  xfi, 

<9B 

(ii)  -7^-  +  V x  E  =  0  in  (0,  T)  xfi, 

(iii)  V  •  D  =  p  in  (0,  T)  x  O, 

(!) 

(iv)  V  •  B  =  0  in  (0,  T)  x  O, 

(v)  E  x  n  =  0  on  (0,  T)  x  3D, 

(vi)  E(0,x)  =  0,  H(0,x)  =  0  in  O. 

The  vector  valued  functions  E  and  H  represent  the  strengths  of  the  electric  and  magnetic  fields, 
respectively,  while  D  and  B  are  the  electric  and  magnetic  flux  densities,  respectively.  The  conduc¬ 
tion  current  density  is  denoted  by  Jc,  while  the  source  current  density  is  given  by  Js.  The  scalar 
p  represents  the  density  of  free  electric  charges  unaccounted  for  in  the  electric  polarization.  We 
assume  perfect  conducting  boundary  conditions  (l,v),  on  the  boundary  3D,  with  unit  outward  nor¬ 
mal  n.  We  also  assume  zero  initial  conditions  for  all  the  unknown  fields.  System  (1)  is  completed 
by  constitutive  laws  that  embody  the  behaviour  of  the  material  in  response  to  the  electromagnetic 
fields.  These  are  given  in  (0,  T)  x  11  in  the  form 

(i)  D(f,x)  =  e0er(x)E(f,x)  +  P(t,x), 

„  (ii)  B(t,x)  =  //0H(t,x). 


(2) 
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For  the  media  that  is  of  interest  to  us,  we  can  neglect  magnetic  effects;  we  also  assume  that 
Ohms’s  law  governs  the  electric  conductivity,  i.e., 


Jc(f,x)  =  <r(x)E(f,x)  in  (0,T)  x  17. 


(3) 


We  will  modify  system  (1)  and  the  constitutive  laws  (2)  by  performing  a  change  of  variables  that 
renders  the  system  in  a  form  that  is  convenient  for  analysis  and  computation.  From  (l,i)  we  have 


—  +  J  Jc(s, x)  ds  )  —  VxH  =  Js  in  (0,  T)  x  17. 


Then  we  define  the  new  variable 


(4) 


D(f,x)  =  D(f,x)  +  /  Jc(s,x)  ds. 


(5) 


Using  definition  (5)  in  (4)  and  (1)  we  obtain  the  modified  system 


(i) 

3D 

~dt 

3B 

(ii) 

~di 

(hi) 

V-l 

(iv) 

V-1 

(v) 

E  x 

(vi) 

E(0. 

(6) 


We  note  that  equation  (6,iii)  follows  from  the  continuity  equation  ^  +  V  •  Jc  =  0,  the  assumption 
that  p(0)  =  0,  and  the  assumption  that  V  •  Js  =  0  (in  the  sense  of  distributions).  The  modified 
constitutive  law  (2,i)  after  substitution  of  (3)  and  (5)  becomes 


D(f, x)  =  eoer(x)E(t, x)  +  /  cr(x)E(s, x)  ds  +  P(7, x). 

Jo 


(7) 


To  describe  the  behaviour  of  the  media’s  macroscopic  electric  polarization  P,  we  employ  a 
general  integral  representation  model  in  which  the  polarization  explicitly  depends  on  the  past 
history  of  the  electric  field.  This  convolution  model  is  sufficiently  general  to  include  microscopic 
polarization  mechanisms  such  as  dipole  or  orientational  polarization  as  well  as  ionic  and  electronic 
polarization  and  other  frequency  dependent  polarization  mechanisms.  The  resulting  constitutive 
law  can  be  given  in  terms  of  a  polarization  or  displacement  susceptibility  kernel  v  as 


P(f,x) 


s ,  x)E(s,  x)  ds. 


(8) 


Thus  the  modified  constitutive  laws  are 


(i)  D(7,  x)  =  e0er(x)E(7,  x)  +  f  {<r(x)  +  u(t  -  s,  x)}  E(s,  x)  ds, 

Jo 

„  (ii)  B(f, x)  =  /i0H(f, x), 


(9) 
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where,  in  the  above  and  henceforth  we  have  dropped  the  ~  symbol  over  D,  at  the  same  time  keeping 
in  mind  that  D  in  definition  (5)  is  the  modified  electric  flux  density.  Let  us  define  the  vector  of 
fields 

u  =  (Ul,u2)T  =  (E,H)T  €  kL1,1(0,  T ;  R6)),  (10) 


and  the  operator 


which  from  (9)  can  be  written  as: 


Lu(i,  x)  = 


D(i,x) 

B  (Lx) 


Lu(t,x)  = 


e0er(x)I3  03  1  /  E(t,  x 


O3  At 0I3  J  \  H(i,x 


t  (\  cj(x)I3  03  1  f  v(t-s,x) I3  0, 


O3  03 


O3  o3 


We  label  the  three  6x6  coefficient  matrices  in  (12)  as 


A(x)  = 


e0er(x)I3  03 
O3  M0I3 


E(s,x) 

H(s,x) 


B(x)  = 


cr(x)I3  03 

O3  03 


c  (t,x)  = 


1 fit,  x)I3  03 

O3  03 


where,  in  the  above  definitions  In  is  an  n  x  n  identity  matrix  and  0n  is  an  n  x  n  matrix  of  zeros, 
n  G  N.  Using  these  definitions  we  may  rewrite  (12)  as 

Lu(t,  x)  =  A(x)u(i,  x)  +  f  B(x)u(s,  x)  ds  +  f  C(i  —  s,  x)u(s,  x)  ds.  (14) 

Jo  Jo 


Next,  we  define  the  Maxwell  operator  M  as 


Mu(i,  x)  =  M 


E(i,x) 

H(t,x) 


VxH(t,x) 
— VxE(t,  x) 


and  the  vector  Js  as 


J  s(t)  =  -Js(t)e  1 


where,  ei  =  (1,  0,  0,  0,  0,  0)T  G  R6,  is  a  unit  basis  vector.  Thus  Maxwell’s  equation  can  be  rewritten 
in  the  form 

d 

(i)  —  Lu  =  Mu  +  Js(i)  in(0,T)xU, 

(ii)  u(0,x)  =  0  in  U,  (1^) 


(iii)  uj  (t,  x)  x  n(x)  =  0  on  (0,  T)  x  9U, 
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where  L  is  the  operator  associated  with  the  constitutive  law  (14),  and  M  is  the  Maxwell  operator 
(15).  Note  that  the  exterior  source  term  Js  has  only  one  non-zero  component. 

We  assume  that  the  structure  that  occupies  the  domain  17  entails  periodic  microstructures 
leading  to  matrices  A,  B  and  C  with  spatially  oscillatory  coefficients.  Specifically,  we  will  assume 
that  er,  a  and  v  are  rapidly  oscillating  spatial  functions. 


Y“  Cell 


a  i 


Figure  1:  Periodic  composite  material  presenting  a  circular  microstructure  with  periodicity  a.  The 
figure  shows  a  decreasing  from  left  to  right. 


3  The  Homogenized  Problem 


The  theory  presented  in  this  section  is  based  on  results  from  [6] .  We  denote  by  Ya  the  reference  cell 
of  the  periodic  structure  that  occupies  17  (see  Figure  1).  The  construction  of  the  homogenized  prob¬ 
lem  involves  solving  for  the  corrector  subterms  £  H^er{Y\  R2),  £  VF1,:L(0,  T;  Hpei(Y]  M2)) 

and  w'(  £  W2,1(0,  T;  iPper(Y;  M2)),  solutions  to  the  corrector  equations 


(i)  /  A(y)Vj,w£  •  Vyv(y)  dy  =  -  /  A(y)efc  •  Vyv(y)  dy, 


>Y 


IY 


(ii)  |  A(y)  Vywfc(i,  y )  +  J  (B(y )  +  C(i  “  y)}  Vyw  k(s,  y)  dsj  •  Vyv(y)dy 
=  -  jv  {B(y)  +  C (t,  y)}  {ek  +  Vyw^}  •  Vyv(y)  dy, 

(iii)  Jy  |  A(y) Vywg(t,  y )  +  J  iB(y)  +  C(f  ~  y)}  Vyw °k(s,  y)ds|  •  Vyv(y)dy 

[  =  -  Jy  A(y)efc  •  Vyv(y)  dy, 


(18) 


for  all  v  £  Hper(Y]M?)  and  k  £  {1,,*-. .  ,6}.  Here,  Hper(Y)  denotes  the  space  of  periodic  functions 
with  vanishing  mean  value.  Note  that  if  we  decompose  v  into  civq  +  C2V2  where  Vi  =  [u  1 , 0]  and 
V2  =  [0,^2],  then  each  set  of  equations  above  decouple  into  an  equation  for  w^i  and  one  for  ~wk,2- 
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For  example,  equation  (18, i)  becomes 


(ia) 

/  An(y )VyWy  •  Vyvi(y)  dy  =  - 

/  An(y)eM  •  Vyvi(y)  dy 

Jy 

r 

Jy 

r 

(19) 

(ib) 

J  A22(y)Vyw(4i2  •  Vyv2(y)  dy  =  - 

j  A22(y)ek}2  •  Vyv2(y)  dy. 

where  An  denotes  the  first  3x3  block  of  A  and  e^i  =  [1, 1, 1,  0,  0,  0]T,  and  similarly  A22  denotes 
the  last  3x3  block  of  A  with  ekt2  =  [0,  0,  0.1, 1, 1]T . 


The  first  corrector  term  u,  from  the  two-scale  expansion  (see  [6]) 

u"  =  u(x)  +  Vyu(x,  y)  +  . . . ,  x  G  11, y  G  Y, 

is  given  as 


(20) 


u(t,x,y)  =  w^4(y)wfc(t,x)  +  /  wk(t  -  s,  y)uk(s,  x)  ds  +  wg(t,  y)'u°(x).  (21) 

Jo 

with  u  G  VF2,1(0,  T;  H^er(Y\  M2)),  where  we  have  considered  the  decompositions  u(t,  x)  =  uk(t,x)ek 
and  u°(x)  =  u9(x)ek,  k  =  1, ...  6.  Rewriting  (21)  in  matrix  form  we  have 

u(t,x,y)  =  wA(y)u(t,x)  +  [  w(f  -  s, y)u(s, x)  ds  +  w°(t, y)u°(x).  (22) 

Jo 

where  w'4  G  M2x6  with  columns  {w^}?=1.  Similarly  w°,w  G  R2x6.  The  expansion  (20)  can  be 
written  as 


=  EX1  +  dwih(x,  y)  +  •  •  • 

Ex2  =  e*2  +  dy2ui(x, y)  +  . . . 

Ex3  =  e*3  +  9y3ni(x,  y)  +  . . . 

HX3  =  Hxi  +  dyiu2(x,  y)  +  •  •  • 

Hx2  =  hx2  +  dy2  fx2(x,  y)  +  •  •  • 

Hx3  =  hx3  +  dy3 U2(x,  y)  + - 

We  now  define  a  new  operator 

£u(t,x)  =du(f,x)  +  8  /  u(s,x)ds+  /  C(t  —  s)u(s,  x)  ds, 

Jo  Jo 


(23a) 

(23b) 

(23c) 

(23d) 

(23e) 

(23f) 

(24) 


where  the  6x6  matrices  A,  B  and  C  are  computed  using  the  solutions  of  system  (18)  as  follows 
(i)  Ak=  [  A(y)  {efc  +  VywA4(y)}dy, 


IY 


(ii)  Bk=  /  B(y)  {ek  +  Vyw(4(y)}dy, 


IY 


(iii)  Ck{t)  =  /  C(*,y){efc  +  Vyw£(y)}dy  +  /  A(y)Vywfc(t,  y)  dy 


(25) 


IY 


IY 


+ 


y  Jo 


{B(y)  +  C (t  -s,y)}  Vyw k(s,  y)  ds  dy, 
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for  k  =  1,2, ...  ,6,  and  where  Ak,Bk,Ck  are  the  fcth  columns  of  the  homogenized  matrices  A,B, 
and  C,  respectively.  In  the  homogenized  problem,  the  electromagnetic  field  u  is  the  solution  of  the 
system 

d 

(i)  — £u  =  Mu  J s  in  (0)  T)  x  P, 

(ii)  u(0,x)  =  0  in  P,  (26) 

(iii)  ui(t)  x  n  =  0  on  (0,  T)  x  9P, 

where  Js  is  as  defined  in  (16),  M  is  as  defined  in  (15),  and  £  is  as  defined  in  (24).  We  note  that 

if  the  initial  conditions  are  nonzero,  then  there  is  a  supplementary  source  term  77°  that  should  be 

introduced  in  the  right  side  of  (26, i),  which  is  given  to  be 

J°(t,x)  =  ug(x)^  ^AVyw k(t)  +  fQ  (B  +  C(f-s))Vyw£(s)  ds^  |  ,  (27) 

k  =  1, . . . ,  6.  See  [6]  for  details. 


4  Reduction  to  Two  Spatial  Dimensions 


We  now  assume  our  problem  possesses  uniformity  in  the  spatial  direction  y  (see  Figure  2  for 
a  schematic  of  the  computational  domain).  In  this  case  Maxwell’s  equations  decouple  into  two 
different  modes,  the  transverse  electric  (TE)  and  transverse  magnetic  (TM)  modes.  Here,  we  are 
interested  in  the  TEy  mode.  The  TEy  mode  involves  the  components  Ex,  Ez  for  the  electric  field 
and  the  component  Hy  of  the  magnetic  field.  Let  x,  y  e  M3  with  x  =  (aq,  aq,  aq)  and  y  =  (yi,  2/2, 2/3)  - 
We  will  use  x  £  R3  for  points  on  the  macro  scale,  and  y  gK3  for  points  on  the  micro  scale  (reference 
cell).  Since  we  assume  uniformity  in  the  aq-direction,  we  may  take  all  derivatives  with  respect  to 
aq  (or  2/2)  to  be  zero.  Then  equation  (17, i)  can  be  written  in  scalar  form  as 


1 

Cl 

1 _ 

1 

H03 

CO 

H* 

to 

1 

_ 1 

dtVX  2 

dx^Hxl  dXiHx  3 

dtDXz 

9x\  HX2 

dtBXl 

— dX3  Ex  2 

dtBX2 

dx^Exi  H-  ^x\EX3 

1 

S’ 

to 

CO 

1 _ 

1 

! 

to 

1 _ 

with  Lu  =  (D,  B)T  Recall  here  that  D  is  the  modified  electric  flux  density,  where  we  have  dropped 
the  ~  notation.  We  may  decouple  system  (28)  into  the  TE  mode, 


1 

s 

C 

dxsHX2  Js 

dtBX2 

= 

dx^Exi  H-  dx\EX3 

ADxs. 

1 

> 

to 

1 _ 

(29) 
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Figure  2:  The  computational  domain  with  surrounding  perfectly  matched  absorbing  layers  (PML). 
The  figure  shows  an  antenna  in  the  region  [x'i, x'2]  x  zc,  X\,X2  6  Xq,  and  the  composite  material 
displaying  periodic  circular  microstructures  in  the  region  Xq  x  Zp.  The  PML  layers  are  terminated 
by  perfectly  conducting  boundary  conditions  (PEC  walls). 


and  the  TM  mode 


5^ 

to 

a 

_ 1 

1 

1 

CO 

tO 

_ 1 

dtDx  2 

= 

&X3  HX1  9X1HX^ 

to 

CO 

1 _ 

1 

1 

to 

1 _ 

(30) 


We  assume  that  our  pulse  is  polarized  to  only  have  an  xi-component.  In  this  case  the  component 
that  is  of  interest  in  our  problem  is  the  EXl  component. 


5  Homogenization  Model  in  Two  Dimensions 


In  a  similar  manner  to  the  three  dimensional  case,  we  may  construct  matrices  ArE,  B1E,  and  C1E 
that  represent  the  constitutive  relations  in  two  dimensions.  Thus  the  constitutive  matrices  are 


(31) 

•  (32) 


Ate  = 

1 

> 

t;i-9 

M 

o' 

;  Bte  = 

r-RTE 

o' 

;  CTE  = 

rnTE 

'“’11 

o' 

L  0 

M  o_ 

0 

0 

0 

0 

A™  (x)  = 


e0er(x)  0 

0  eoer(x) 


B™(X)  = 


<r(x)  0 

0  <r(x) 


Cu(t,x)  = 


Ktx) 

0 


0 

i/(t,x) 


The  homogenized  solution  for  the  TE  mode  is  obtained  from  the  formal  asymptotic  expansion 
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(23)  as 


EXl  =  EXl  +  dyi ui  (x,  y)  +  . . .  (33a) 

Ex3  =  ex3  +  dy3u i  (x,  y )  +  . . .  (33b) 

Hx2  =  HX2  +  dy2u2(x,  y)  + -  (33c) 

Also,  since  we  are  assuming  uniformity  in  the  x2  direction,  we  set  the  term  dy2u2(x.,y )  to  zero.  So 
(33)  becomes 


EX1  =  Ex i  +  dyiui(x,  y)  +  . . .  (34a) 

EX3  =  Ex3  +  dy3ui  (x,  y )  +  . . .  (34b) 

HZ2  =  HX2.  (34c) 

Hence  the  homogenized  electric  field  for  the  TE  mode  is 

EQ  =  E  +  Vy«i(x,y)  +  . . . ,  (35) 


where  the  gradient  operator  in  this  case  is  Vy  =  (dyi ,  dy3)T .  Therefore  we  only  need  to  solve  for 
tti(x,  y),  which  in  turn  only  depends  on  the  first  component  of  w^4,  w)1 ,  and  w^,  for  k  =  1,2;  we 
will  refer  to  these  as  ,  wy ,  and  wk,  respectively. 


Let  us  again  denote  by  Y  the  reference  cell  of  the  periodic  structure  that  occupies  H  C  M2. 
The  construction  of  the  two-dimensional  homogenized  problem  involves  solving  for  the  corrector 
subterms  £  ifp0r(T ;  M),  wk  £  1T1,1(0,  T;  H^er{Y ;  M))  and  w £  £  1T2,1(0,  T;  H^er.{Y\  M)),  solutions 
to  the  corrector  equations 


(i)  /  AiT(y)vywfc  •  Vyu(y)  dy  =  -  /  A1i1E(y)efc  •  Vyu(y)  dy, 


'Y 


IY 


(ii)  [  Aj^(y)Vywk(t,y)  ■  Vyu(y)dy 
JY 


+ 


Y  JO 


{B11  (y)  +  Cn(t  -  S,y)}  Vyivk(s,  y)  ds  •  Vyu(y)dy 


=  -  /  {Bn  (y)  +  cn  (^y)}{efc  +  VyWfc} -Vyu(y)  dy, 


IY 


(iii)  J  A-(y)Vy  ,o(,y)  .  Vyu(y)dy 


+  /  /  {BnE(y)  +  CnE(t-s,y)}  Vyt^(s,y)ds- Vyu(y)dy 

IY  Jo 

=  ~  j  A41b(y)efc  •  Vyu(y)  dy, 


(36) 


for  all  v  £  Hpei(Y;M),  k  =  1,2  and  ei  =  [l,0]T,e2  =  [0, 1]T.  Since  the  initial  conditions  that  we 
have  chosen  are  zero,  we  will  not  need  to  calculate  the  corrector  w9.  We  will  only  need  to  solve  for 
and  wk.  Once  we  have  solved  for  these  corrector  terms,  we  can  then  construct  the  homogenized 


11 


matrices  from 

(i)  (^nE)fc  =  Jy  A^E(y)  {efc  +  Vyt()^(y)}  dy, 

(ii)  (Bj (E)fe  =  J ^  BjIE(y )  {efc  +  Vyu>£(y)}  dy, 

(iii)  (Cj^)k(t)  =  [  CjE(t,y)  {ek  +  Vyrc^(y)}  dy  +  [  A^E(y)Vyri)fc(t, y)  dy 


(37) 


+  /  /  {Bn(y)  +  C^E(f  -  s,y)}  Vywk(s,y)  ds  dy, 

Jy  Jo 

where  e^,  k  =  1,  2  are  the  basis  vectors  in  M2,  (A™)*,,  (7?™)*,,  and  (C™)*,  are  the  fcth  columns  of 
the  matrices  AjE,Bji‘,  and  CjE,  respectively,  and  the  homogenized  matrices  are  given  by 


(38) 


ate  = 

-A™ 

O' 

;  Bte  = 

OS 

M 

oi 

;  CTE  = 

'cJF  O' 

0 

1 

O 

O 

1 _ 

O 

O 

_ 1 

The  corresponding  system  of  equations  in  the  TE  mode  are 

d 


(i)  — £1Ev  =  Mltjv  +  Jg  in  (0,  T )  x  17, 

(ii)  v(0,x)  =  0  in  17, 

l  (iii)  v3nxi  -  v  177.3,3  =  0  on  (0,  T)  x  <917, 


(39) 


where  v  =  {EXl ,  HX2 ,  EX3)T,  n  =  (nXl ,  nX3)T  is  the  unit  outward  normal  vector  to  <917,  the  operator 
£te  is  defined  as 


£teu(7,x)  =  ATEu(f,x)  +  B1E  I  u(s,x)  ds  +  f  CTE(t  —  s)u(s,  x)  ds, 

Jo  Jo 


(40) 


and  M1E  is  the  two-dimensional  curl  operator 


0  -<9 


X3 


a 


X3 


0  -a 


X\  I  ! 


(41) 


0  a 


Xl 


and  jJE(f)  =  —  Jsei,  ei  =  [1,0]T. 

6  Models  for  Polarization 


The  constitutive  law  in  (8)  is  sufficiently  general  to  include  models  based  on  differential  equations 
and  systems  of  differential  equations  or  delay  differential  equations  whose  solutions  can  be  ex¬ 
pressed  through  fundamental  solutions  (in  general  variation-of-paranreters  representations)-see  [2] 
for  details.  A  number  of  known  polarization  laws  can  be  readily  treated. 


1.  The  choice  of  the  kernel  function 


zz(t,x) 


e0(es(x)  -  eoo(x)) 
t(x) 


-t/r(x) 

7 


x  £  Xq  X  Zd 


(42) 
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in  the  dielectric  Xq  x  Zp  (see  Figure  2)  corresponds  to  the  differential  equation  of  the  Debye 
Model  for  orientational  or  dipolar  polarization  given  by 

tP*  +  P  =  e0(es  -  eoo)E,  (43) 

D(t,x)  =  e0eoo(t,  x)E(f,  x)  +  P(t,x).  (44) 


Here  es  is  the  static  relative  permittivity.  The  presence  of  instantaneous  polarization  is  ac¬ 
counted  for  in  this  case  by  the  coefficient  in  the  electric  flux  equation.  That  is,  er  =  e^ 
in  Xq  x  Zd,  er  =  1  in  air.  The  remainder  of  the  electric  polarization  is  seen  to  be  a  decaying 
exponential  with  relaxation  parameter  r,  driven  by  the  electric  field,  less  the  part  included 
in  the  instantaneous  polarization.  This  model  was  first  proposed  by  Debye  [11]  to  model  the 
behaviour  of  materials  that  possess  permanent  dipole  moments.  The  magnitude  of  the  polar¬ 
ization  term  P  represents  the  degree  of  alignment  of  these  individual  moments.  The  choice  of 
coefficients  in  (43)  gives  a  physical  interpretation  to  es  and  as  the  relative  permittivities  of 
the  medium  in  the  limit  of  the  static  field  and  very  high  frequencies,  respectively.  In  the  static 
case,  we  have  Pt  =  0,  so  that  P  =  eo(es  —  Coo)E  and  D  =  eseoE.  For  very  high  frequencies, 
rPt  dominates  P  so  that  P  ~  0  and  D  =  CoofoE. 

2.  The  general  model  also  includes  the  Lorentz  model  for  electronic  polarization  which,  in  dif¬ 
ferential  form,  is  represented  with  the  second  order  equation: 

P«H — Pt  +  WqP  =  eoWpE,  (45) 

T 

D  =  eocCoE  +  P .  (46) 

In  (45),  Up  is  called  the  plasma  frequency  and  is  defined  to  be 

uP  =  u0y/es  -  eoo,  (47) 


where  uq  is  the  resonance  frequency.  A  simple  variation  of  constants  solution  yields  the 
correct  kernel  function 

v{t)  =  °  pe~f//2rsin  (48) 

D) 


vq  = 


(49) 


3.  For  more  complex  dielectric  materials,  a  simple  Debye  or  Lorentz  polarization  model  is  often 
not  adequate  to  characterize  the  dispersive  behaviour  of  the  material.  One  can  then  turn 
to  combinations  of  Debye,  Lorentz,  or  even  more  general  nth  order  mechanisms  as  well  as 
Cole-Cole  type  (fractional  order  derivatives)  models.  We  again  refer  the  reader  to  [2,  10]  for 
details. 
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7  Numerical  Discretization  for  the  Cell  Problem 

7.1  Spatial  discretization  via  finite  elements 

We  divide  the  reference  cell  Y  into  elementary  rectangles,  and  consider  7/,  to  be  a  uniform  mesh 
with  elements  {K}  of  edge  length  h.  We  define  the  finite  dimensional  space 

=  {vh  |  vh  G  C°(Y),vh\K  G  Qi  for  all  K  G  7^},  (50) 

which  approximates  171(Y).  In  (50),  the  space  Q\  is  defined  as  Q i  =  Pn,  where,  for  £q,  k2  G  NU{0} 

Tfcife  =  {p(xi,x2)b(xi,x2)  =  ^  ^  a,ijx\x32,  a,ij  G  M}.  (51) 

0<i<fci  0<j<k2 

Thus,  p  i  is  the  space  of  the  bilinear  functions  of  two  variables,  and  V/,  is  the  space  of  continuous 
piecewise  bilinear  functions.  We  now  consider  the  subspace 

Vper,h  =  {Vh  |  Vh  G  Vh,Vh(yi,0)  =  Vh{y i,l)  and  vh(0 ,y2)  =  vh(l ,y2)  for  all  yi,y2  G  [0,1]},  (52) 

of  Vh,  where  Y  =  [0, 1]  x  [0, 1]. 

Thus  the  spatially  discrete  problem  is  to  find  wkh  G  Vperi/j,  u>k,h.  e  W1,1(0,  T;  Vper,/i)  and 
w^h  ^  IY2,1(0, T;  Vper^),  solutions  to  the  corrector  equations 

(i)  e0  j r  er,h(y) •  VyVh{y)  dy  =  -e0  J ^  er>h(y)ek  ■  VyVh(y)  dy, 

(h)  eo  J  er^h(y)VyWkth(t,y)  -VyVh(y)dy 

+  J  J  Wh{y)  +  ish{t- s,y)}VyWk,h{s,y) -VyVh{y)  ds  dy 

=  -  J^Wh(y)  +  vh(t,y)}  {efc  +  vjtwjfft}  -VyVh(y)  dy,  (53) 

(hi)  eo  J  er>h(y)V$wl'h(t,y)  •  VyD/l(y)dy 

+  J  J  {<7h(y)  +  Mt  -  s,  y)>  VyW^s,  y )  •  vjufc(y)  ds  dy 

=  — e0  j  er,h{y)ek  ■  VyVh(y)  dy, 

for  all  Vh  G  Vperh  and  k  =  1,2,  and  ei  =  [l,0]T,e2  =  [0, 1]T,  and  Vy  is  a  discrete  approximation 
to  the  gradient.  In  (53),  we  have  used  the  definitions  (32)  of  the  matrices  A™,  B™  and  C™. 


7.2  Time  discretization  via  a  recursive  convolution  approach 

Since  the  susceptibility  kernel  u(t,  x),  is  exponential  in  nature  for  many  materials  of  interest,  we 
can  use  recursion  to  compute  the  discretized  time  convolution  of  the  susceptibility  kernel  with  the 
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electric  field  in  the  corrector  subproblems,  e.g.,  (53, ii)  and  (53, hi).  A  similar  approach,  known  as 
the  recursive  convolution  (RC)  method  has  been  used  to  compute  the  discrete  convolution  terms 
that  appear  in  Maxwell’s  equations  [18,  15].  Let  tn  =  nAt,  for  some  timestep  At  and  let  Vn  denote 
the  time  component  V(nAt)  for  any  vector  field  V.  Thus,  the  equation  (53, ii)  can  be  discretized 
in  time  as 


eo  /  er,/i(y)V£itf£fc(y)  •  VyVh(y)dy 


IY 


+ 


y 

nAt 


Y  JO 


lY 


Wh( y)  +  Vh{nAt  -  s,  y)}  Vy wk,h(s,  y)  •  VyVh(y)  ds  dy 
W(y)  +  Vh(nAt,  y)}  |efc  +  v£u)(Jh}  '  ^yVh{y)  dy. 


(54) 


We  will  assume  that  all  field  components  are  constant  over  each  time  interval  of  length  At.  Thus, 
assuming  that  all  field  vectors  are  zero  for  t  <  0  we  have 


eo  /  er,h(y)v}w2h(y)  ■  Vjuft(y)dy 


IY 


+ 


n—  1 

E 


K(y)  +  vh(nAt  -s,y)}  Vjwg^y)  •  V$vh(y)  dsdy 


JY  m= 0  JmAt 

=  ~  J  {o-fc(y)  +  ^(nAt,  y)}  |efc  +  Vy-w^h]  ■  Vyhh(y)  dy. 
Interchanging  the  space  and  time  integral  terms  and  rearranging  we  obtain 
eo  er,h( y)V}u}%h(y)  •  Vjvh(y)dy 


n  1  p  (  /*(m+l)At  \ 

+  E  /  /  - 5- y> ds  vZ<k'(y>  ■  vy^(y>  dy 

m= 0Jy  VmAt  J 

n—1  „ 

+  /,  Hy)  •  Vyh/l(y)  dy 

rn=0  '  Y 

W(y)  +  ^.(nAt,y)}  jefc  +  VyW^h]  •  VyVh{y)  dy. 


Let  us  define 


*C(y)  = 

Using  definition  (57)  in  (56)  we  have 


r»(m+l)  At 


'  m  At 


vh{nAt  -  s,y)  ds. 


71—1 


eo  er,h{ y)Vjw£)fc(y)  •  Vjufc(y)dy  +  ^  /  CV^JV)  ■  v£uh(y)  dy 

*^y  m=0JY 

n—1  .. 

+  ^2At  ah{y)^y™k£l(y)  ■  v^h(y)  dy 

m=0 

{cr/,(y)  +  i//»(nAt,y)}  |efc  +  •  V$vh(y)  dy. 


(55) 


(56) 


(57) 


(58) 
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Upon  rearranging,  we  have 

eo  j r  {er,h{ y)  +  A t  crh( y)  +  v™_  1}  V}w%h(y)  ■  Vfch(y)dy 

_ 2  TL 2 

+  E  /  <E<E(y)  •  dy  +  E  At  j ,  ^(y)E<f(y)  •  v$«fc(y)  dy  (59) 

m=0  'A  ?n=0  'A 

=  -  Wh{ y)  +  Vh{nAt,  y)}  |efc  +  Vju)jjh}  •  VyVh(y)  dy. 

In  equation  (59)  all  the  terms  can  be  computed  with  the  knowledge  of  just  the  solution  at  time 
tn  =  nAt,  i.e. ,  Vy w^h,  except  the  discrete  convolution  term  involving  for  m  =  0, ...  ,n  —  2. 
We  now  show  that  the  discrete  convolution  of  all  previous  Vy field  values  and  the  discrete 
susceptibility  function  can  be  reduced  to  recursive  updating  of  a  single  vector  on  each  element  in 
the  finite  element  mesh,  which  involves  a  matrix  vector  multiplication  at  each  time  step. 


7.3  Recursive  convolution  for  Debye  polarization 


We  consider  the  case  of  Debye  polarization  in  the  remainder  of  this  paper.  In  this  case  the  spatially 
discrete  susceptibility  function  is  defined  as 


vh{t,y)  =  eo(e^(y)  e°°,ft(y))e-t/T*(y),  yeY,t>0, 

w) 

and  from  (57)  and  (60)  the  function  is  defined  as 

i/£  (y)  =  [ {m+1)At  eo(es,h{y)  ~  eooM)  e-(tn-s)/TH(y)  dg>  y  G  Y. 

Jm.At  W) 

From  (59)  we  define  the  summation 

n—  2  „ 

n  =  E  /  <^y<hl(y)  •  Vj«ft(y)  dy,  K  G  Th, 

m=0  JK 


(60) 


(61) 


(62) 


where  n”  is  defined  in  (61).  Equation  (62)  constitutes  the  time  discrete  convolution  of  the  suscep¬ 
tibility  function  Uh  and  all  field  values  of  Vy w™h  up  to  the  nth  time  step. 

Then  for  the  n  +  1  st  step  we  find 

71—1 


n+1  =  E  /  "m'VXdl  y)  •  v}®»(y)  dy 

771=0  ^  K 


71—2 


(63) 


E  /  .  <+1Vjffi”+1(y)  •  Vjsj(y)  dy  +  Ctjvjffi^fy)  •  Vjs,,(y)  dy. 
_ n  J  K  J  K 


711=0 


From  (61)  we  can  derive  the  identity 

<+1(y)  =  ^(y)«-A,,Thiy}- 


(64) 
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Using  (64)  in  (63)  we  obtain 


n— 2 


t?1  =  E  /  .  •  Vjsk(y)  dy 

m=0*'A 


(65) 


+  C+|vXJy)'V^(y)  dy. 
J  K 


In  the  spatial  discretization  of  the  reference  cell  Y",  we  will  assume  that  the  functions  spatially 
dependent  t> t(y),  eSj^(y),  eoo,/i(y)>  and  <J/l( y)  are  constant  on  each  element  K  of  the  finite  element 
mesh  7^.  With  this  assumption,  for  every  element  77  £  7),.  we  have 


n— 2  „ 

=  g—At/Th(K)  J2  /  ^  i&Vjwg^y)  •  Vjuh(y)  dy 

m=0  ^  K 


(66) 


+  ClVj^’V^y)  dy. 

J  K 

Combining  (66)  and  (62),  we  obtain  the  recursion 

n+1  =  e-^^n  +  [  K-Wy<ht y)  •  vy^(y)  dy’ 

J  K 

where  n™ji((77),  from  the  definition  (61),  can  be  calculated  to  be 

^(tf)  =  eo (ea>h(K)  -  e^h{K))e^2nAt^K\eAt/T^  -  1). 


(67) 


(68) 


The  finite  element  function  'IT  which  is  defined  to  be  on  element  K  in  the  triangulation  7), 
can  be  calculated  as 


n—1 


K+'  =  E  / .  C+1  Vjr0”+‘(y)  .  Vjs*(y)  dy, 

m=0JY 


n—1  r. 

=  E  E  /  <+1Vj»^+1(y)  ■  V}sft(y)  dy,  (69) 

A'erh  m=o  ^  x 

=  E  (e-A,/™<*'ln+  /  CIUUW y)-Vjsk(y)dy). 

xerh  A  Jh  ‘ 

Thus  we  have  the  recursion 

=  E  +  [  ^Vj<>A(y)  •  VyUft(y)  dy.  (70) 

A"eTh 

We  note  that  given  the  value  4/)^  for  every  element  K  e  7^  at  tn  =  nAt,  we  can  calculate  'h)(+1  by 
equation  (70)  which  only  requires  the  knowledge  of  the  solution  vector  w^h  at  the  nth  step!  The 
solutions  wTh  for  m  <  n  are  not  needed  to  be  stored  in  memory  for  this  calculation.  We  can  treat 
the  discrete  convolution  in  (53,iii)  in  the  same  manner  as  above. 
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The  calculation  of  the  discrete  homogenized  matrix  also  requires  the  computation  of  a 
time  integral  which  at  time  tn  =  nAt  requires  the  knowledge  of  the  solutions  w™h  at  previous 
times  rriAt.  m  =  1, . .  .  n.  Thus,  we  again  try  to  obtain  a  recursive  formula  that  will  aid  in  the 
computation  of  these  time  integrals.  Once  we  have  solved  for  the  discrete  corrector  terms,  we  can 
then  construct  the  discrete  homogenized  matrices  from 

(  (i)  MnEft)fe  =  Aii,h(y)  {efc  +  y)}  dy, 

(ii)  (BTi Eh)k  =  [  Bjg,(y)  lek  +  Vjuijjh(y)}  dy, 

Jy  (71) 

(iii)  {CjEh)k{t)  =  Cjg^y)  {efc  +  VjtDj£fc(y)}  dy  +  Aj^h{y)V$wk,h{t,y)  dy 


+  /  /  {Bn ,h(y)  +  cn ,h(t  -  a,y)}  vju?fc,h(s,y)  ds  dy, 


Y  JO 


where,  as  in  (37),  ek,  k  =  1,2  are  the  basis  vectors  in  M2,  (AjEh)k,  (BjEh)k,  and  {CjEh)k  are 
the  kth  columns  of  the  discretized  matrices  AjEh,  l3jEh,  an(i  CjEh,  respectively,  and  the  discrete 
homogenized  matrices  are  calculated  as 


AlE  = 


AjEh  O' 
0  no 


t3lE  = 


Bu Eh  o' 
0  0 


CjE  = 


1 


o' 

0  0 


(72) 


Using  the  definitions  of  the  discrete  matrices  we  rewrite  (71)  as 

(i)  MnEft)fe  =  eo  eoo,h( y)  {efc  +  Vjwjfft(y)}  dy, 

(ii)  =  [  <rh(y)  {efc  +  Vju>j£h(y)}  dy, 


(iii)  (CiT5Jfc(y  =  J  vh{t,y)  {efc  +  Vjtuj£fc(y)}  dy  +  e0^eOO)fe(y)V^fc,/l(i,y)  dy 


(73) 


+ 


nAt 


Y  JO 


W(y)  +  Vh(nAt  -  s,  y)}  Vy wkjh(s,  y)  ds  dy. 


The  term 


nAt 


Vh{nAt  -  s,y)VyWkh(s,y)  ds  dy, 


(74) 


Y  JO 


in  equation  (73, iii)  can  be  computed  using  a  recursive  formula  that  does  not  require  knowledge  of 
the  solutions  w™h  at  times  mAt ,  m  =  1, . . . ,  n.  To  derive  a  recursive  formula,  which  again  need  to 
use  properties  of  the  susceptibility  function  .  We  define  the  term  in  (74)  to  be 


71—1 


TJi  =  E  /  ^y<h  W)dy. 

m= 0Jy 


(75) 


where  is  defined  in  (61).  We  then  have 

n  « 

r,;+1  =  E  /  <.+1vJ»S1(».y)  dy 

m-0Jy 


71—1 


Y  JYU™+1V iWS'1(S’y^  dy  +  JYUn+ly/y'WkX1(s,y)  dy 


m= 0 


(76) 
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Using  (64)  in  (76)  we  obtain  the  recursion 


rpn- (-1  _ 

1h  ~ 


71—1 

E 

m= 0 


g-At/rhiy)^ 


'Y 


dy  +  X^v^C.y)  dy 


=  E  e-MlT{K)mK)  +  /  ^+1V^<+1(S,y)dy, 

AGT,  ^ 


where 


n— 1 

T£(A1  =  E  /  <V^-+1(S,y)  dy. 
rn=0JK 


(77) 


(78) 


8  Numerical  Examples 


In  this  section  we  present  results  of  initial  computations  based  on  the  above  discretization  schemes 
that  we  have  carried  out  for  both  circular  and  square  microstructures.  In  these  computations 
we  have  used  a  conjugate  gradient  method  to  solve  the  resulting  linear  systems  that  arise  after 
discretizing  our  model  in  space  and  time. 


8.1  Examples  with  varying  relative  permittivity 


In  the  following  examples  we  choose  the  value  of  the  relative  permittivity  to  be 


et  =  1.003,  if  x  E  S, 
ee  =  2.7,  if  x  €  Y/S, 

where  S  is  the  microstructure  that  is  enclosed  inside  the  reference  cell.  These  values  are  taken  from 
experimental  measurements  for  air  and  polyurethane  material  which  are  the  primary  components 
of  the  insulating  foam  described  in  the  introduction. 


8.1.1  Circular  microstructures 


We  first  consider  the  example  of  a  composite  material  which  possesses  circular  microstructures  in 
two-dimensions.  We  will  solve  a  cell  problem  in  the  reference  cell  Y  =  [0, 1]  x  [0, 1],  in  which  the 
relative  permittivity  is  given  in  equation  (79)  where  S  is  the  circular  microstructure  enclosed  inside 
the  reference  cell,  as  depicted  in  Figure  3.  In  this  test  case  we  will  assume  that  es  =  e^,  a  =  0  and 
r  =  3.16  x  10~8.  Hence  the  parameters  a  and  r  are  constant  over  the  entire  dielectric  material.  Since 
eoo  =  es j  we  have  that  v{t)  =  0  for  all  time  t.  Thus  in  this  example,  the  model  possesses  instant 
polarization  but  does  not  have  a  hysteretic  term  in  the  polarization.  Our  numerical  simulation  is 
performed  on  a  51  x  51  nodes  mesh  grid.  We  define  the  inclusion  volume  fraction  /  as  the  ratio 

area  of  inclusion 

j  =  - — - —  =  area  of  inclusion 

area  of  domain  r 


(80) 
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Y“  Cell 


Figure  3:  (left)  Periodic  composite  material  presenting  a  circular  mi¬ 
crostructure  with  periodicity  a.  (right)  The  reference  cell  Y  =  [0, 1]  x 
[0,1]- 


For  /  =  0.5  the  homogenized  matrix  *4™  is 

1.68009  -2.92  x  10“4 

-2.92  x  10“4  1.68009 


aTE  -  c 


and  the  homogenized  matrices  Bjf'  and  CjfJ  are  the  zero  matrices. 


(81) 


In  the  numerical  solution  for  the  homogenized  model  we  have  approximated  the  circular  rni- 
crostructure  in  a  staircase  fashion.  In  Table  1  we  present  the  homogenized  relative  permittivities 
for  different  inclusion  volume  fractions  /,  and  different  refinements  of  the  mesh  grid  that  is  imposed 
on  the  reference  cell  Y .  Here  h  denotes  the  mesh  step  size.  In  this  table  we  note  that  the  lower 
and  upper  off-diagonal  entries,  (A i^hod;  and  (»4ii£')uodj  respectively,  of  the  homogenized  matrix 
decrease  at  least  as  fast  as  O(h)  as  the  mesh  grid  is  refined.  The  error  in  the  off-diagonal 
terms  is  probably  due  to  the  inaccurate  representation  of  the  circular  microstructure  by  a  staircase 
approximation.  Hence  Aj^  is  approximately  diagonal  with  identical  diagonal  entries.  For  the  case 
of  /  =  0.5  we  have 


Ajp  «  1.68  eo 


1.0  0 
0  1.0 


(82) 


Thus  the  homogenized  or  effective  relative  permittivity  for  /  =  0.5  is  er  =  1.68. 


In  Figure  4  we  plot  the  relative  effective  permittivity  versus  the  inclusion  volume  fraction  for 
the  method  discussed  in  this  paper  and  compare  it  to  different  theoretical  mixture  models,  namely 
the  Maxwell-Garnett  formula,  the  Bruggeman  mixture  rule,  as  well  as  the  weighted  average  of  the 
different  relative  permittivities.  (We  refer  the  reader  to  [22]  for  a  discussion  of  the  various  theoretical 
mixing  formulas.)  We  note  how  close  our  predicted  values  of  the  effective  relative  permittivity  are 
to  those  predicted  by  these  mixing  rules.  The  prediction  of  the  effective  relative  permittivity  of  the 
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/ 

/sc 

h 

(^)d 

("^ll^Olod 

(^l/')uod 

0.7854 

0.78 

0.02 

1.23068 

5.92955  x  10”3 

5.92955  x  10“3 

0.78 

0.01 

1.23186 

1.86948  x  10-3 

1.86948  x  10-3 

0.785 

0.005 

1.23347 

1.16484  x  10~3 

1.16484  x  10~3 

0.785 

0.0025 

1.23767 

5.41779  x  10-4 

5.41779  x  10-4 

0.5027 

0.502 

0.02 

1.68009 

-2.92023  x  10-4 

-2.92023  x  10-4 

0.502 

0.01 

1.68299 

-7.31937  x  10“5 

-7.31937  x  10“5 

0.502 

0.005 

1.68509 

-1.94965  x  10-5 

-1.94965  x  10-5 

0.502 

0.0025 

1.68523 

-5.02621  x  10“6 

-5.02621  x  10“6 

0.7124 

0.7124 

0.02 

1.34434 

-8.69546  x  10~19 

-2.42287  x  10-17 

0.7129 

0.01 

1.34789 

-3.21656  x  10-18 

-1.73589  x  10-17 

0.7124 

0.005 

1.35081 

2.19197  x  10“16 

5.59808  x  10"16 

0.7125 

0.0025 

1.35216 

-1.96101  x  10“18 

5.536  x  10”16 

0.1257 

0.1252 

0.02 

2.40606 

-6.64092  x  10-4 

-6.64092  x  10-4 

0.1253 

0.01 

2.40637 

-1.72749  x  10“4 

-1.72749  x  10“4 

0.1255 

0.005 

2.40609 

-4.39996  x  10“5 

-4.39996  x  10“5 

0.1255 

0.0025 

2.40626 

-1.15057  x  10-5 

-1.15057  x  10-5 

Table  1:  The  effective  relative  permittivities  for  different  volume  frac¬ 
tions  with  a  circular  microstructure.  The  table  lists  the  diagonal  entries, 
(ylf]£')d,  of  the  homogenized  matrix  as  well  as  the  lower  and  upper 
off-diagonal  entries  and  (v4|",E)uod:  respectively,  for  different 

levels  of  refinement  of  the  mesh  grid  on  the  reference  cell  Y .  Here  /  is 
the  volume  fraction  and  fsc  is  the  computed  volume  fraction  using  the 
staircase  approximation. 
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Figure  4:  Effective  relative  permittivity  distribution  for  random  mixtures  with 
circular  microstructures. 


composite  mixture  eeg  according  to  these  formulas  is  given  by 

eeg  =  ee  +  2/ee - - — — - -  (Maxwell-Garnett),  (83) 

ei  +  ee-  J{ei  -  €e) 

(1  _  /) €e  6eff  +  fe‘  e°fi  =  0  (Bruggeman).  (84) 

€e  T"  6eff  "T"  ^eff 

These  formulas  hold  for  randomly  distributed  circular  inclusions  of  permittivity  in  a  homogeneous 
environment  of  permittivity  ee.  The  inclusions  occupy  a  volume  fraction  /  of  the  homogeneous 
medium.  These  different  mixing  models  predict  different  effective  permittivity  values  for  a  given 
mixture.  There  are  also  bounds  that  limit  the  range  of  the  predictions.  These  bounds  are  the 
Weiner  bounds  given  by 


eeff.max  =  /ej  +  (1  -  /)ee  (weighted  average), 


(85) 


^eff,min 


f€e  +  (1  -  fW 


(86) 


These  bounds  hold  for  all  values  of  e,;  and  ee.  Other  mixing  models  such  as  power-law  models,  the 
Lichtenecker  formula,  etc.,  also  exist  in  the  literature.  We  again  refer  the  reader  to  [23,  22]  for 
further  details.  In  Figure  5  we  plot  the  solution  vectors  Wi  and  for  /  =  0.5,  over  the  reference 
cell  Y .  In  Figure  6  we  plot  a  topview  of  the  corresponding  solution  vectors  in  the  domain  Y. 
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0.08, 


Figure  5:  Plot  of  the  solution  vectors  Wi  (top)  and  W2  (bottom)  for  an  inclu¬ 
sion  fraction  volume  /  =  0.5  on  a  50  x  50  cells  mesh  grid.  The  homogenized 
value  of  er  for  this  case  is  ~  1.68 


axis 
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Figure  6:  Top  view  of  the  solution  vectors  wf  (top)  and  (bottom) 
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8.1.2  Square  microstructures 


We  next  replaced  the  circular  microstructures  in  two-dimensions  with  square  microstructures  and 
repeated  the  experiment  performed  in  Section  8.1.1.  We  solved  a  cell  problem  in  the  reference  cell 
Y  =  [0, 1]  x  [0, 1]  in  which  the  value  of  the  relative  permittivity  is  given  in  (79)  where  S  is  now 
a  square  microstructure  that  is  enclosed  inside  the  reference  cell.  In  this  test  case,  as  before,  we 
will  assume  that  es  =  Coo,  c  =  0  and  r  =  3.16  x  10~8.  Thus,  the  parameters  a  and  r  are  constant 
over  the  entire  dielectric  material.  Since  =  es,  we  again  have  that  v(t)  =  0  for  all  time  t.  Our 


Y“  Cell 


Figure  7:  (left)  Periodic  composite  material  presenting  a  square  mi¬ 
crostructure  with  periodicity  a.  (right)  The  reference  cell  Y  =  [0, 1]  x 
[0,1]- 


numerical  simulation  is  performed  on  a  51  x  51  nodes  mesh  grid.  For  / 
matrix  *4™  is 


A 


TE 

11 


eo 


1.65057  -3.04533  x  10~18 

-1.8576  x  10'18  1.65057 


0.52  the  homogenized 


(87) 


and  the  homogenized  matrices  and  CjE  are  the  zero  matrices.  In  this  case  the  sides  of  the 
square  are  aligned  with  the  mesh  grid  on  the  reference  cell  Y ,  so  we  do  not  have  the  staircase  error 
that  was  observed  in  the  case  of  circular  microstructures.  Again  AjE  is  approximately  diagonal 
with  identical  diagonal  entries.  For  the  case  of  /  =  0.52  we  have 


Aii  «  1.65  e0 


1.0  0 
0  1.0 


(88) 


Thus  the  homogenized  or  effective  relative  permittivity  for  /  =  0.52  is  er  =  1.65.  In  Figure  8  we  plot 
the  effective  relative  permittivity  for  different  volume  fractions  and  compare  these  values  to  those 
obtained  for  circular  microstructures.  We  note  here  how  close  the  effective  relative  permittivity 
values  are  for  circular  and  square  microstructures  of  the  same  volume  fraction.  In  Figure  9  we  plot 
the  solution  vectors  wf  and  w for  /  =  0.52,  over  the  reference  cell  Y .  In  Figure  10  we  plot  a 
topview  of  the  corresponding  solution  vectors  in  the  domain  Y . 
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8.2  Example  with  varying  relaxation  times 
8.2.1  Circular  microstructures 


We  now  consider  an  example  of  a  composite  material  which  possesses  circular  microstructures  in 
which  the  value  of  the  relaxation  time  r  is  given  to  be 


r(x) 


Tj  =  1.58  x  10  8,  if  x  €  S, 

Te  =  3.16  x  10— 8 ,  if  x  G  Y/S. 


(89) 


In  this  test  case,  we  will  assume  that  the  other  parameters  namely,  es,  a  and  are  constant  over 
the  entire  dielectric  material  with  Coo  =  5.5,  es  =  78.2  and  a  =  1.0  x  10~5. 


In  Figure  11  we  plot  the  results  of  the  cell  problem.  This  figure  plots  the  effective  dielectric 
response  function  (DRF)  as  a  function  of  time.  We  compare  the  homogenized  DRF  with  four 
different  cases.  The  low  frequency  case  corresponds  to  one  in  which  the  relaxation  time  r  is  given 
to  be  the  weighted  average  of  r*  and  re 

Tow  =  fn  +  (1  -  f)re.  (90) 


Figure  8:  Comparison  of  the  effective  relative  permittivity  distribution  for 
square  and  circular  microstructures 
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y  axis 


0.08. 


y  axis 


Figure  9:  Plot  of  the  solution  vectors  wf  (top)  and  (bottom)  for  a  square 
inclusion  fraction  volume  /  =  0.52  on  a  50  x  50  cells  mesh  grid.  The  homoge¬ 
nized  value  of  er  for  this  case  is  ~  1.65 


axis 
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Figure  10:  Top  view  of  the  solution  vectors  wf  (top)  and  w £  (bottom)  for 
square  microstructures. 
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The  high  frequency  case  corresponds  to  the  relaxation  time 


Thigh 


1  -/ 
Te 


-1 


(91) 


The  low  and  high  frequency  averages  for  the  relaxation  time  r  were  observed  in  numerical  experi¬ 
ments  that  were  performed  in  [3] ,  in  which  a  probabilistic  approach  is  taken  to  an  electromagnetic 
interrogation  problem  of  a  dielectric  medium  that  is  a  mixture  of  two  Debye  media  with  different 
relaxation  times  r.  The  low  frequency  average  corresponds  to  the  average  relaxation  time  that 
is  observed  when  the  composite  media  is  interrogated  by  electromagnetic  pulses  at  frequencies  of 
106  Hz,  whereas  the  high  frequency  average  was  observed  at  higher  frequencies  of  109  —  1011  Hz 
(frequencies  are  considered  low  here  if  the  corresponding  angular  frequency  lo  <  1  jr).  In  Figure  11 
we  also  compare  the  homogenized  DRF  to  the  weighted  average  of  the  corresponding  DRF’s,  i.e., 

^wa(i)  =  M  +  (1  -  f)Vei  (92) 


where 


(i)  Vi(t) 

(ii)  ve{t) 


_  £o(es  eOO )  —t/Ti 

Ti 

_  €o(es  —  eoo)  — t/Te 

6  j 

Te 


(93) 


as  well  as  to  another  average  that  we  call  the  inverse  weighted  average,  which  is  the  inverse  of  the 
weighted  averages  of  the  inverses  of  the  two  DRF’s,  i.e., 


(1  -/) 

Ve 


-1 


(94) 


where  v,  and  ve  are  defined  in  (93, ii),  and  (93,iii),  respectively.  Clearly  Figure  11  demonstrates 
good  agreement  between  the  homogenized  DRF  and  the  other  approximations.  Finally  in  Figure 
12  we  plot  the  homogenized  DRF  as  a  function  of  the  inclusion  volume  fraction  at  times  t  =  0.556 
ns  and  t  =  27.78  ns.  For  t  «  min(Tj,Te)  one  can  show  that  viWSL(t)  ~  i '(t;r\ow)  as  is  clearly  seen 
in  the  top  plot  of  this  figure.  Similarly  we  have  vwa(t)  ~  v(t;  Thigh)-  We  observe  that  for  all  values 
of  /  the  plot  of  the  homogenized  DRF  lies  between  the  other  approximations  for  both  small  and 
large  times. 


8.2.2  Square  microstructures 


We  next  repeated  the  above  numerical  calculation  for  square  microstructures  in  which  the  value  of 
the  relaxation  time  r  is  given  to  be 


r(x) 


Ti  =  1.58  x  10  8,  if  x  e  S, 

Te  =  3.16  x  10  8 ,  if  x  G  Y/S. 


(95) 


As  before  we  have  assumed  that  the  other  parameters  namely,  es,  a  and  are  constant  over  the 
entire  dielectric  material  with  =  5.5,  es  =  78.2  and  cr  =  1.0  x  10~5.  In  Figure  13  we  plot  the 
homogenized  DRF  as  a  function  of  the  inclusion  volume  fraction  at  time  t  =  0.556  ns  (top)  and 
time  t  =  27.78  ns  (bottom).  In  Figure  14  we  plot  the  homogenized  DRF  as  a  function  of  time,  for 
/  =  0.52.  In  both  Figures  13  and  14  we  observe  that  the  difference  between  the  values  of  the  DRF 
for  square  versus  circular  microstructures  is  of  the  order  of  10~4  —  10~5. 
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Figure  11:  (top)  The  effective  dielectric  response  function  (DRF)  u(t)  as  a 
function  of  time,  for  volume  fraction  /  =  0.5.  We  compare  the  homogenized 
DRF  with  the  DRF  that  is  a  weighted  average  of  the  corresponding  DRF’s,  as 
well  averages  that  are  observed  in  other  experiments,  (bottom)  Semilog  plot 
of  the  DRF’s. 
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Figure  12:  (top)  The  effective  dielectric  response  function  (DRF)  v(t)  as  a 
function  of  the  inclusion  volume  fraction  for  t  =  0.556  ns  and  (bottom)  t  = 
27.78  ns. 
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Figure  13:  (top)  The  effective  dielectric  response  function  (DRF)  v{t)  in  the 
case  of  square  microstructures  for  three  different  values  of  the  volume  fraction 
/.  (bottom)  The  difference  in  the  effective  DRF  between  square  and  circular 
microstructures.  In  all  three  cases  the  differences  are  of  the  order  of  10-5,  while 
we  note  from  the  top  plot  that  the  functions  themselves  are  of  magnitude  10~2. 
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Figure  14:  (top)  The  effective  dielectric  response  function  (DRF)  u(t)  as  a 
function  of  the  inclusion  volume  fraction  /  for  t  =  0.556  ns  and  (bottom) 
t  =  27.78  ns.  The  difference  between  the  values  of  the  DRF  for  the  square 
and  the  circular  microstructures  are  on  the  order  of  10~4  —  10-5. 
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9  Conclusions  and  Future  Directions 

Taking  into  account  the  large  deformations  of  the  corrector  functions  (see  Figures  5,  6,  9,  and 
10)  we  can  expect  that  the  correction  of  the  electric  and  magnetic  fields  as  given  by  (20)  will  be 
of  crucial  interest.  The  results  of  the  cell  problem  with  circular  microstructures  for  calculating 
effective  relative  permittivities  agree  well  with  the  mixing  formulae  results  that  are  available  in  the 
literature.  However,  there  are  several  advantages  to  taking  a  homogenization  approach  to  obtain 
effective  dielectric  parameters.  First,  mixing  formulae  are  not  available  for  general  microstructure 
geometries.  Close  examination  of  the  materials  of  long  term  interest  to  us  reveal  that  the  foam 
consists  of  hexagonal-like  cells  in  a  microstructure  configuration  that  is  95%  gas  surrounded  by 
a  matrix  of  polyurethane  which  has  an  estimated  relative  permittivity  of  2.7.  A  homogenization 
approach  will  permit  us  to  ascertain  the  sensitivity  of  the  effective  permittivities  with  respect  to  the 
micro-geometry.  While  we  have  used  a  staircase  method  to  approximate  the  circular  microstructure 
in  the  cell  problem,  the  approach  also  permits  better  approximations  (e.g.,  a  fictitious  domain 
approach)  for  the  circular  (and  other  shape)  microstructures  in  the  cell  problem. 

A  second  advantage  of  a  homogenization  approach  is  the  flexibility  it  affords  in  assumptions 
about  material  polarization  laws.  Here  in  proof-of-concept  calculations  we  have  used  a  Debye 
medium  for  polarization  but  could  with  equal  ease  investigate  the  cases  of  Lorentz  media,  compos¬ 
ite  laws,  and  other  higher  order  dispersive  media.  Moreover,  as  we  investigate  other  mechanisms, 
it  may  be  computationally  advantageous  to  treat  directly  the  polarization  law  as  a  side  constraint. 
That  is,  the  hysteretic  term  of  the  polarization  is  most  generally  represented  by  an  integral  term 
which  involves  a  dielectric  response  function.  Our  numerical  approximation  to  this  approach  in¬ 
volved  a  recursive  convolution  method  to  approximate  the  integral  term.  Our  current  efforts  entail 
an  approach  in  which  the  homogenized  model  involves  an  ordinary  differential  equation  for  the 
hysteretic  part  of  the  polarization  term.  This  approach  should  facilitate  numerical  approximation 
of  the  solution  of  the  electromagnetic  interrogation  problem  for  a  wide  class  of  assumed  polarization 
mechanisms. 

The  initial  results  of  the  cell  problems  for  the  calculation  of  the  dielectric  response  function 
compare  well  with  the  results  that  were  obtained  in  [3],  where  a  probabilistic  approach  is  taken  to 
an  electromagnetic  interrogation  problem  in  a  dielectric  medium  that  is  a  mixture  of  two  Debye 
media  with  different  relaxation  times.  Each  approach  has  its  own  conceptual,  theoretical  and 
computational  advantages  that  merit  further  comparisons. 
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